Note: Open scripts.Rproj first, then script. To easily use relative paths, click the down button next to knit and then click “Knit Directory –> Project Directory”. This should make loading and saving files much easier.
This script is based off of Langfelder P, Horvath S (2008) WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 2008, 9:559 (link to paper) Red text indicates regions that require the user to modify. Regardless, the user should check over all code blocks to ensure that everything is running correctly.
Load packages
library(tidyverse, quietly = TRUE)
library(genefilter, quietly = TRUE)
library(edgeR, quietly = TRUE)
library(RColorBrewer, quietly = TRUE)
library(WGCNA, quietly = TRUE)
library(flashClust, quietly = TRUE)
library(gridExtra, quietly = TRUE)
library(ComplexHeatmap, quietly = TRUE)
library(goseq, quietly = TRUE)
library(dplyr, quietly = TRUE)
library(clusterProfiler, quietly = TRUE)
library(simplifyEnrichment, quietly = TRUE)
Change file names and conditions where appropriate.
# Input unfiltered data
treatmentinfo.file <- "samples.tsv"
gcount.file <- "salmon.numreads.tsv"
# Output DEG results
WGCNA_results.file <- "salmon.numreads.WGCNA_results.tsv"
WGCNA_samples.file <- "salmon.numreads.WGCNA_samples.tsv"
# pOverA gene filtering thresholds (only process genes with over these cutoffs)
pOverA.cutoff.P <- 0.8 # >80% of samples
pOverA.cutoff.A <- 10 # >10 read counts
# Column name to use for analysis of data using WGCNA
group <- "group"
# Minimum module size - I chose 30 as it is the default value chosen by most studies using WGCNA.
minModuleSize <- 30
Load the input file containing the treatment information.
treatmentinfo <- read.csv(treatmentinfo.file, header = TRUE, sep = "\t", fileEncoding="UTF-8-BOM") # Read in file
head(treatmentinfo)
## sample_id unit lib_type fq1 fq2 species plug_id
## 1 Pacuta_ATAC_TP1_1043 1 pe SRR14610932 SRR14610932 Pacuta 1043
## 2 Pacuta_ATAC_TP1_1775 1 pe SRR14610931 SRR14610931 Pacuta 1775
## 3 Pacuta_ATAC_TP1_2363 1 pe SRR14610930 SRR14610930 Pacuta 2363
## 4 Pacuta_ATAC_TP3_1041 1 pe SRR14610929 SRR14610929 Pacuta 1041
## 5 Pacuta_ATAC_TP3_1471 1 pe SRR14610928 SRR14610928 Pacuta 1471
## 6 Pacuta_ATAC_TP3_1637 1 pe SRR14610927 SRR14610927 Pacuta 1637
## treatment timepoint reef ploidy group
## 1 ATAC 1 Lilipuna.Fringe 3 Group4
## 2 ATAC 1 Reef.42.43 2 Group5
## 3 ATAC 1 Reef.18 3 Group3
## 4 ATAC 3 Reef.11.13 3 Group2
## 5 ATAC 3 Reef.35.36 3 Group3
## 6 ATAC 3 Reef.11.13 3 Group2
# Check we have the right column names
headers <- c("sample_id")
if( all(headers %in% colnames(treatmentinfo)) ){
cat(paste(treatmentinfo.file, "has the required columns:", paste(headers, collapse=', ')))
} else {
stop(paste(treatmentinfo.file, "is missing required columns:", paste(headers, collapse=', ')))
}
## samples.tsv has the required columns: sample_id
Filter treatmentinfo to just the samples that we want to
analyze. Generally we would use all samples for WGCNA, however, in some
cases we might want to split the data if we have samples from different
experiments. Adjust filtering as
required.
#treatmentinfo <- treatmentinfo %>%
# filter(condition %in% c("ATAC", "HTAC") & timepoint=="8") %>%
# mutate(condition = factor(condition, levels = c("ATAC", "HTAC")))
head(treatmentinfo)
## sample_id unit lib_type fq1 fq2 species plug_id
## 1 Pacuta_ATAC_TP1_1043 1 pe SRR14610932 SRR14610932 Pacuta 1043
## 2 Pacuta_ATAC_TP1_1775 1 pe SRR14610931 SRR14610931 Pacuta 1775
## 3 Pacuta_ATAC_TP1_2363 1 pe SRR14610930 SRR14610930 Pacuta 2363
## 4 Pacuta_ATAC_TP3_1041 1 pe SRR14610929 SRR14610929 Pacuta 1041
## 5 Pacuta_ATAC_TP3_1471 1 pe SRR14610928 SRR14610928 Pacuta 1471
## 6 Pacuta_ATAC_TP3_1637 1 pe SRR14610927 SRR14610927 Pacuta 1637
## treatment timepoint reef ploidy group
## 1 ATAC 1 Lilipuna.Fringe 3 Group4
## 2 ATAC 1 Reef.42.43 2 Group5
## 3 ATAC 1 Reef.18 3 Group3
## 4 ATAC 3 Reef.11.13 3 Group2
## 5 ATAC 3 Reef.35.36 3 Group3
## 6 ATAC 3 Reef.11.13 3 Group2
Load the input file containing the gene count matrix.
gcount <- as.data.frame(read.csv(gcount.file, header = TRUE, sep = "\t", fileEncoding="UTF-8-BOM")) # Read in file
# Check we have the right column names
headers <- c("Name", treatmentinfo$sample_id)
if( all(headers %in% colnames(gcount)) ){
cat(paste(gcount.file, "has the required columns:", paste(headers, collapse=', ')))
} else {
stop(paste(gcount.file, "is missing required columns:", paste(headers, collapse=', ')))
}
## salmon.numreads.tsv has the required columns: Name, Pacuta_ATAC_TP1_1043, Pacuta_ATAC_TP1_1775, Pacuta_ATAC_TP1_2363, Pacuta_ATAC_TP3_1041, Pacuta_ATAC_TP3_1471, Pacuta_ATAC_TP3_1637, Pacuta_ATAC_TP4_1060, Pacuta_ATAC_TP4_1762, Pacuta_ATAC_TP4_2002, Pacuta_ATAC_TP5_1059, Pacuta_ATAC_TP5_1563, Pacuta_ATAC_TP5_1757, Pacuta_ATAC_TP6_1050, Pacuta_ATAC_TP6_1468, Pacuta_ATAC_TP6_1542, Pacuta_ATAC_TP7_1047, Pacuta_ATAC_TP7_1445, Pacuta_ATAC_TP7_2413, Pacuta_ATAC_TP8_1051, Pacuta_ATAC_TP8_1755, Pacuta_ATAC_TP8_2012, Pacuta_ATAC_TP9_1141, Pacuta_ATAC_TP9_1594, Pacuta_ATAC_TP9_2357, Pacuta_ATAC_TP10_1159, Pacuta_ATAC_TP10_1559, Pacuta_ATAC_TP10_1641, Pacuta_ATAC_TP11_1103, Pacuta_ATAC_TP11_1777, Pacuta_ATAC_TP11_2306, Pacuta_ATHC_TP1_1207, Pacuta_ATHC_TP1_2743, Pacuta_ATHC_TP1_2977, Pacuta_ATHC_TP3_1219, Pacuta_ATHC_TP3_2534, Pacuta_ATHC_TP3_2750, Pacuta_ATHC_TP4_1220, Pacuta_ATHC_TP4_2733, Pacuta_ATHC_TP4_2993, Pacuta_ATHC_TP5_1296, Pacuta_ATHC_TP5_2212, Pacuta_ATHC_TP5_2877, Pacuta_ATHC_TP6_1254, Pacuta_ATHC_TP6_2870, Pacuta_ATHC_TP6_2999, Pacuta_ATHC_TP7_1281, Pacuta_ATHC_TP7_2409, Pacuta_ATHC_TP7_2878, Pacuta_ATHC_TP8_1459, Pacuta_ATHC_TP8_2564, Pacuta_ATHC_TP8_2861, Pacuta_ATHC_TP9_1451, Pacuta_ATHC_TP9_2873, Pacuta_ATHC_TP9_2979, Pacuta_ATHC_TP10_1205, Pacuta_ATHC_TP10_2197, Pacuta_ATHC_TP10_2550, Pacuta_ATHC_TP11_1147, Pacuta_ATHC_TP11_2668, Pacuta_ATHC_TP11_2879, Pacuta_HTAC_TP1_1653, Pacuta_HTAC_TP1_2005, Pacuta_HTAC_TP1_2414, Pacuta_HTAC_TP3_1617, Pacuta_HTAC_TP3_1642, Pacuta_HTAC_TP3_2026, Pacuta_HTAC_TP4_1581, Pacuta_HTAC_TP4_1701, Pacuta_HTAC_TP4_1767, Pacuta_HTAC_TP5_1303, Pacuta_HTAC_TP5_1571, Pacuta_HTAC_TP5_1707, Pacuta_HTAC_TP6_1330, Pacuta_HTAC_TP6_1466, Pacuta_HTAC_TP6_1744, Pacuta_HTAC_TP7_1487, Pacuta_HTAC_TP7_1728, Pacuta_HTAC_TP7_2072, Pacuta_HTAC_TP8_1329, Pacuta_HTAC_TP8_1765, Pacuta_HTAC_TP8_2513, Pacuta_HTAC_TP9_1302, Pacuta_HTAC_TP9_1486, Pacuta_HTAC_TP9_1696, Pacuta_HTAC_TP10_1225, Pacuta_HTAC_TP10_1536, Pacuta_HTAC_TP10_2064, Pacuta_HTAC_TP11_1582, Pacuta_HTAC_TP11_1596, Pacuta_HTAC_TP11_1647, Pacuta_HTHC_TP1_1239, Pacuta_HTHC_TP1_1676, Pacuta_HTHC_TP1_2210, Pacuta_HTHC_TP3_1227, Pacuta_HTHC_TP3_1418, Pacuta_HTHC_TP3_2527, Pacuta_HTHC_TP4_1169, Pacuta_HTHC_TP4_1343, Pacuta_HTHC_TP4_2195, Pacuta_HTHC_TP5_1168, Pacuta_HTHC_TP5_1415, Pacuta_HTHC_TP5_2087, Pacuta_HTHC_TP6_1138, Pacuta_HTHC_TP6_1595, Pacuta_HTHC_TP6_1721, Pacuta_HTHC_TP7_1090, Pacuta_HTHC_TP7_1427, Pacuta_HTHC_TP7_1820, Pacuta_HTHC_TP8_1184, Pacuta_HTHC_TP8_1709, Pacuta_HTHC_TP8_2304, Pacuta_HTHC_TP9_1131, Pacuta_HTHC_TP9_2202, Pacuta_HTHC_TP9_2305, Pacuta_HTHC_TP10_1238, Pacuta_HTHC_TP10_1732, Pacuta_HTHC_TP10_2300, Pacuta_HTHC_TP11_1416, Pacuta_HTHC_TP11_2185
# Cleanup gcount data
gcount <- gcount %>%
column_to_rownames("Name") %>% # Makes "Name" column rownames
round() %>% # Round
select(treatmentinfo$sample_id) # Select just columns in filtered treatmentinfo file
# View dataset attributes
d <- dim(gcount); print(paste("rows:",d[[1]]," columns:",d[[2]], sep=''))
## [1] "rows:33259 columns:119"
head(gcount)[,1:3]
## Pacuta_ATAC_TP1_1043
## Pocillopora_acuta_KBHIv2___RNAseq.g16115.t1 78
## Pocillopora_acuta_KBHIv2___TS.g16868.t1 216
## Pocillopora_acuta_KBHIv2___RNAseq.g25384.t1 4
## Pocillopora_acuta_KBHIv2___RNAseq.g30122.t1 0
## Pocillopora_acuta_KBHIv2___RNAseq.g12817.t1 1
## Pocillopora_acuta_KBHIv2___TS.g9097.t1 90
## Pacuta_ATAC_TP1_1775
## Pocillopora_acuta_KBHIv2___RNAseq.g16115.t1 47
## Pocillopora_acuta_KBHIv2___TS.g16868.t1 121
## Pocillopora_acuta_KBHIv2___RNAseq.g25384.t1 0
## Pocillopora_acuta_KBHIv2___RNAseq.g30122.t1 0
## Pocillopora_acuta_KBHIv2___RNAseq.g12817.t1 0
## Pocillopora_acuta_KBHIv2___TS.g9097.t1 12
## Pacuta_ATAC_TP1_2363
## Pocillopora_acuta_KBHIv2___RNAseq.g16115.t1 32
## Pocillopora_acuta_KBHIv2___TS.g16868.t1 211
## Pocillopora_acuta_KBHIv2___RNAseq.g25384.t1 3
## Pocillopora_acuta_KBHIv2___RNAseq.g30122.t1 0
## Pocillopora_acuta_KBHIv2___RNAseq.g12817.t1 2
## Pocillopora_acuta_KBHIv2___TS.g9097.t1 107
Filter gcount to just genes with enough data across
columns (exclude genes which are mostly zeros)
# Create filter for the counts data
filt <- filterfun(pOverA(pOverA.cutoff.P, pOverA.cutoff.A))
gfilt <- genefilter(gcount, filt)
# Identify genes to keep by count filter
keep <- gcount[gfilt,]
# Identify gene lists
keep <- rownames(keep)
# Gene count data filtered in PoverA, P percent of the samples have counts over A
gcount_filt <- as.data.frame(gcount[which(rownames(gcount) %in% keep),])
d <- dim(gcount); print(paste("gcount - rows:",d[[1]]," columns:",d[[2]], sep='')) # Before filtering
## [1] "gcount - rows:33259 columns:119"
d <- dim(gcount_filt); print(paste("gcount_filt - rows:",d[[1]]," columns:",d[[2]], sep='')) # After filtering
## [1] "gcount_filt - rows:15467 columns:119"
Determine library size.
libSize.df <- data.frame(libSize=colSums(gcount_filt))
Make DGE object.
DGEdat <- DGEList(counts=as.matrix(gcount_filt),
samples=treatmentinfo,
group=treatmentinfo[[group]])
dim(DGEdat$counts)
## [1] 15467 119
DGEdat <- calcNormFactors(DGEdat)
head(DGEdat$samples)
## group lib.size norm.factors sample_id unit
## Pacuta_ATAC_TP1_1043 Group4 4128082 0.9769418 Pacuta_ATAC_TP1_1043 1
## Pacuta_ATAC_TP1_1775 Group5 4557568 0.9699660 Pacuta_ATAC_TP1_1775 1
## Pacuta_ATAC_TP1_2363 Group3 4629399 0.9728287 Pacuta_ATAC_TP1_2363 1
## Pacuta_ATAC_TP3_1041 Group2 4870724 0.9988819 Pacuta_ATAC_TP3_1041 1
## Pacuta_ATAC_TP3_1471 Group3 3699133 0.9679178 Pacuta_ATAC_TP3_1471 1
## Pacuta_ATAC_TP3_1637 Group2 3459979 0.9964851 Pacuta_ATAC_TP3_1637 1
## lib_type fq1 fq2 species plug_id treatment
## Pacuta_ATAC_TP1_1043 pe SRR14610932 SRR14610932 Pacuta 1043 ATAC
## Pacuta_ATAC_TP1_1775 pe SRR14610931 SRR14610931 Pacuta 1775 ATAC
## Pacuta_ATAC_TP1_2363 pe SRR14610930 SRR14610930 Pacuta 2363 ATAC
## Pacuta_ATAC_TP3_1041 pe SRR14610929 SRR14610929 Pacuta 1041 ATAC
## Pacuta_ATAC_TP3_1471 pe SRR14610928 SRR14610928 Pacuta 1471 ATAC
## Pacuta_ATAC_TP3_1637 pe SRR14610927 SRR14610927 Pacuta 1637 ATAC
## timepoint reef ploidy group.1
## Pacuta_ATAC_TP1_1043 1 Lilipuna.Fringe 3 Group4
## Pacuta_ATAC_TP1_1775 1 Reef.42.43 2 Group5
## Pacuta_ATAC_TP1_2363 1 Reef.18 3 Group3
## Pacuta_ATAC_TP3_1041 3 Reef.11.13 3 Group2
## Pacuta_ATAC_TP3_1471 3 Reef.35.36 3 Group3
## Pacuta_ATAC_TP3_1637 3 Reef.11.13 3 Group2
Log transform the counts matrix for the next plots
DGEdat.cpm <- DGEdat # Make a copy the edgeR dataset
DGEdat.cpm$counts <- cpm(DGEdat.cpm$counts, log=TRUE, prior.count=5) # Log transform the copy for the next plots
gsampleDists_dev <- dist(t(DGEdat.cpm$counts)) # Calculate distance matrix
gsampleDistMatrix_dev <- as.matrix(gsampleDists_dev) # Distance matrix
rownames(gsampleDistMatrix_dev) <- colnames(DGEdat.cpm) # Assign row names
colnames(gsampleDistMatrix_dev) <- NULL # Assign col names
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255) # Assign colors
pheatmap(gsampleDistMatrix_dev, # Plot matrix
clustering_distance_rows=gsampleDists_dev, # Cluster rows
clustering_distance_cols=gsampleDists_dev, # Cluster columns
col=colors) # Set colors
pca <- prcomp(t(DGEdat.cpm$counts)) # Calculate eigengenes
pc.data <- summary(out<-prcomp(t(DGEdat.cpm$counts))); pc.data
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 27.0984 24.7536 23.00132 20.99766 16.05678 14.56444
## Proportion of Variance 0.1363 0.1137 0.09818 0.08182 0.04785 0.03937
## Cumulative Proportion 0.1363 0.2500 0.34817 0.43000 0.47784 0.51721
## PC7 PC8 PC9 PC10 PC11 PC12
## Standard deviation 13.69539 12.15218 11.81621 10.84512 10.4334 10.37330
## Proportion of Variance 0.03481 0.02741 0.02591 0.02183 0.0202 0.01997
## Cumulative Proportion 0.55202 0.57942 0.60534 0.62716 0.6474 0.66733
## PC13 PC14 PC15 PC16 PC17 PC18 PC19
## Standard deviation 9.56842 9.03839 8.90666 8.68261 8.11781 7.84050 7.1939
## Proportion of Variance 0.01699 0.01516 0.01472 0.01399 0.01223 0.01141 0.0096
## Cumulative Proportion 0.68433 0.69949 0.71421 0.72820 0.74043 0.75184 0.7614
## PC20 PC21 PC22 PC23 PC24 PC25 PC26
## Standard deviation 7.03702 6.63725 6.43582 6.07434 6.0078 5.89786 5.6402
## Proportion of Variance 0.00919 0.00818 0.00769 0.00685 0.0067 0.00646 0.0059
## Cumulative Proportion 0.77063 0.77881 0.78649 0.79334 0.8000 0.80649 0.8124
## PC27 PC28 PC29 PC30 PC31 PC32 PC33
## Standard deviation 5.47977 5.41242 5.32948 5.14454 5.00973 4.82876 4.6441
## Proportion of Variance 0.00557 0.00544 0.00527 0.00491 0.00466 0.00433 0.0040
## Cumulative Proportion 0.81797 0.82341 0.82868 0.83359 0.83825 0.84257 0.8466
## PC34 PC35 PC36 PC37 PC38 PC39 PC40
## Standard deviation 4.58795 4.49855 4.38008 4.3454 4.22217 4.13238 4.10966
## Proportion of Variance 0.00391 0.00376 0.00356 0.0035 0.00331 0.00317 0.00313
## Cumulative Proportion 0.85048 0.85424 0.85780 0.8613 0.86461 0.86778 0.87092
## PC41 PC42 PC43 PC44 PC45 PC46 PC47
## Standard deviation 4.08065 4.05213 3.90462 3.83247 3.79697 3.73699 3.72268
## Proportion of Variance 0.00309 0.00305 0.00283 0.00273 0.00268 0.00259 0.00257
## Cumulative Proportion 0.87401 0.87705 0.87988 0.88261 0.88528 0.88788 0.89045
## PC48 PC49 PC50 PC51 PC52 PC53 PC54
## Standard deviation 3.64532 3.62479 3.56088 3.54082 3.51107 3.44885 3.42440
## Proportion of Variance 0.00247 0.00244 0.00235 0.00233 0.00229 0.00221 0.00218
## Cumulative Proportion 0.89291 0.89535 0.89771 0.90003 0.90232 0.90453 0.90670
## PC55 PC56 PC57 PC58 PC59 PC60 PC61
## Standard deviation 3.37977 3.37074 3.33918 3.31446 3.26671 3.24470 3.22975
## Proportion of Variance 0.00212 0.00211 0.00207 0.00204 0.00198 0.00195 0.00194
## Cumulative Proportion 0.90882 0.91093 0.91300 0.91504 0.91702 0.91897 0.92091
## PC62 PC63 PC64 PC65 PC66 PC67 PC68
## Standard deviation 3.22092 3.2022 3.17807 3.16082 3.13362 3.1168 3.09767
## Proportion of Variance 0.00193 0.0019 0.00187 0.00185 0.00182 0.0018 0.00178
## Cumulative Proportion 0.92284 0.9247 0.92661 0.92847 0.93029 0.9321 0.93387
## PC69 PC70 PC71 PC72 PC73 PC74 PC75
## Standard deviation 3.07047 3.04402 3.0250 3.00161 2.98065 2.97374 2.96028
## Proportion of Variance 0.00175 0.00172 0.0017 0.00167 0.00165 0.00164 0.00163
## Cumulative Proportion 0.93562 0.93734 0.9390 0.94071 0.94236 0.94400 0.94563
## PC76 PC77 PC78 PC79 PC80 PC81 PC82
## Standard deviation 2.95028 2.9346 2.92463 2.89990 2.89350 2.86431 2.83641
## Proportion of Variance 0.00162 0.0016 0.00159 0.00156 0.00155 0.00152 0.00149
## Cumulative Proportion 0.94724 0.9488 0.95043 0.95199 0.95354 0.95507 0.95656
## PC83 PC84 PC85 PC86 PC87 PC88 PC89
## Standard deviation 2.81509 2.80595 2.79664 2.76909 2.76666 2.75655 2.73046
## Proportion of Variance 0.00147 0.00146 0.00145 0.00142 0.00142 0.00141 0.00138
## Cumulative Proportion 0.95803 0.95949 0.96094 0.96237 0.96379 0.96520 0.96658
## PC90 PC91 PC92 PC93 PC94 PC95 PC96
## Standard deviation 2.71741 2.70707 2.69606 2.69002 2.67369 2.6500 2.63861
## Proportion of Variance 0.00137 0.00136 0.00135 0.00134 0.00133 0.0013 0.00129
## Cumulative Proportion 0.96795 0.96931 0.97066 0.97200 0.97333 0.9746 0.97592
## PC97 PC98 PC99 PC100 PC101 PC102 PC103
## Standard deviation 2.63269 2.61371 2.60259 2.58884 2.57400 2.5481 2.52325
## Proportion of Variance 0.00129 0.00127 0.00126 0.00124 0.00123 0.0012 0.00118
## Cumulative Proportion 0.97721 0.97848 0.97974 0.98098 0.98221 0.9834 0.98459
## PC104 PC105 PC106 PC107 PC108 PC109 PC110
## Standard deviation 2.50312 2.49367 2.46185 2.44319 2.41778 2.40138 2.38468
## Proportion of Variance 0.00116 0.00115 0.00112 0.00111 0.00108 0.00107 0.00106
## Cumulative Proportion 0.98576 0.98691 0.98804 0.98914 0.99023 0.99130 0.99235
## PC111 PC112 PC113 PC114 PC115 PC116 PC117
## Standard deviation 2.36155 2.33573 2.32859 2.28930 2.26562 2.22176 2.21337
## Proportion of Variance 0.00103 0.00101 0.00101 0.00097 0.00095 0.00092 0.00091
## Cumulative Proportion 0.99339 0.99440 0.99541 0.99638 0.99733 0.99825 0.99916
## PC118 PC119
## Standard deviation 2.12894 1.91e-13
## Proportion of Variance 0.00084 0.00e+00
## Cumulative Proportion 1.00000 1.00e+00
plot(out)
#DGEdat.cpm$samples$timepoint <- gsub("TP", "", DGEdat.cpm$samples$timepoint)
DGEdat_PCcor <- lapply(DGEdat.cpm$samples, as.factor)
#DGEdat_PCcor <- as.tibble(lapply(DGEdat_PCcor, as.numeric))
DGEdat_PCcor <- as.tibble(DGEdat_PCcor)
## Warning: `as.tibble()` was deprecated in tibble 2.0.0.
## ℹ Please use `as_tibble()` instead.
## ℹ The signature and semantics have changed, see `?as_tibble`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
DGEdat_PCcor <- cbind(DGEdat.cpm$samples$lib.size, DGEdat_PCcor)
rownames(DGEdat_PCcor) <- DGEdat.cpm$samples$sample_id
# Make a dataframe containing all plotting info
percentVar <- pca$sdev^2/sum(pca$sdev^2) # Save % variation by PC1 and PC2
d <- data.frame(PC1 = pca$x[, 1],
PC2 = pca$x[, 2],
Group = DGEdat.cpm$samples[[group]],
DGEdat.cpm$samples, name = colnames(DGEdat.cpm$counts)
)
allgenes_PCA <- ggplot(data = d, aes_string(x = "PC1", y = "PC2", colour = "Group")) + # NOTE: can add 'shape' to be appropriate metadata column to add extra info
geom_point(size = 3) +
xlab(paste0("PC1: ", round(percentVar[1] * 100), "% variance")) +
ylab(paste0("PC2: ", round(percentVar[2] * 100), "% variance")) +
coord_fixed() +
theme_bw() + # Set background color
theme(panel.border = element_blank(), # Set border
panel.grid.major = element_blank(), # Set major gridlines
panel.grid.minor = element_blank(), # Set minor gridlines
axis.line = element_line(colour = "black", size = 0.6), # Set axes color
plot.background=element_blank(), # Set the plot background
axis.title = element_text(size = 14), # Axis title size
axis.text = element_blank()) # Axis text size and view plot
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
allgenes_PCA
Transpose the filtered gene count matrix so that the gene IDs are rows and the sample IDs are columns.
datExpr <- as.data.frame(t(DGEdat.cpm$counts))
Check for genes and samples with too many missing values with goodSamplesGenes. There shouldn’t be any because we performed pre-filtering
gsg = goodSamplesGenes(datExpr, verbose = 3)
## Flagging genes and samples with too many missing values...
## ..step 1
gsg$allOK # Should return TRUE if not, the R chunk below will take care of flagged data
## [1] TRUE
Remove flagged samples if the allOK is FALSE
print(paste("Number genes before filtering: ", ncol(datExpr), sep='')) # Number genes before
## [1] "Number genes before filtering: 15467"
if (!gsg$allOK) # If the allOK is FALSE...
{
# Optionally, print the gene and sample names that are flagged:
if (sum(!gsg$goodGenes)>0) {
printFlush(paste("Removing genes:", paste(names(datExpr)[!gsg$goodGenes], collapse = ", ")))
}
if (sum(!gsg$goodSamples)>0) {
printFlush(paste("Removing samples:", paste(rownames(datExpr)[!gsg$goodSamples], collapse = ", ")))
}
# Remove the offending genes and samples from the data:
datExpr = datExpr[gsg$goodSamples, gsg$goodGenes]
}
print(paste("Number genes after filtering: ", ncol(datExpr), sep='')) # Number genes after
## [1] "Number genes after filtering: 15467"
Look for outliers by examining the following dendrogram.
sampleTree = hclust(dist(datExpr), method = "average");
# Plot the sample tree: Open a graphic output window of size 12 by 8 inches
# The user should change the dimensions if the window is too large or too small.
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5, cex.axis = 1.5, cex.main = 2)
Exclude selected samples. We will exlcude Pacuta_HTHC_TP11_1416 and Pacuta_HTAC_TP9_1302 because they appear to be outliers on the dendrogram.
print(paste("Number samples before outlier removal: ", nrow(datExpr), sep='')) # Number samples before
## [1] "Number samples before outlier removal: 119"
row_names_df_to_remove<-c("Pacuta_HTHC_TP11_1416", "Pacuta_HTAC_TP9_1302") # Picked from above plot
datExpr <- (datExpr[!(row.names(datExpr) %in% row_names_df_to_remove),])
treatmentinfo <- (treatmentinfo[!(treatmentinfo$sample_id %in% row_names_df_to_remove),])
print(paste("Number samples after outlier removal: ", nrow(datExpr), sep='')) # Number samples after
## [1] "Number samples after outlier removal: 117"
The soft thresholding power (β) is the number to which the co-expression similarity is raised to calculate adjacency. The function pickSoftThreshold performs a network topology analysis. The user chooses a set of candidate powers, however the default parameters are suitable values.
# Choose a set of soft-thresholding powers
powers <- c(seq(from = 1, to=10, by=0.5)) # Create a string of numbers from 1 through 10, increasing my 0.5 incraments
# Call the network topology analysis function
sft <-pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)
## pickSoftThreshold: will use block size 2892.
## pickSoftThreshold: calculating connectivity for given powers...
## ..working on genes 1 through 2892 of 15467
## Warning: executing %dopar% sequentially: no parallel backend registered
## ..working on genes 2893 through 5784 of 15467
## ..working on genes 5785 through 8676 of 15467
## ..working on genes 8677 through 11568 of 15467
## ..working on genes 11569 through 14460 of 15467
## ..working on genes 14461 through 15467 of 15467
## Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
## 1 1.0 0.0345 0.574 0.970 2820.00 2750.00 4700
## 2 1.5 0.0626 -0.534 0.950 1460.00 1390.00 3100
## 3 2.0 0.3460 -1.110 0.951 817.00 747.00 2170
## 4 2.5 0.6020 -1.460 0.964 486.00 423.00 1600
## 5 3.0 0.7330 -1.710 0.962 304.00 249.00 1210
## 6 3.5 0.8100 -1.830 0.969 198.00 152.00 944
## 7 4.0 0.8440 -1.900 0.966 134.00 95.30 750
## 8 4.5 0.8610 -1.940 0.966 92.70 61.00 605
## 9 5.0 0.8820 -1.940 0.971 66.10 40.00 495
## 10 5.5 0.8910 -1.930 0.974 48.10 26.60 409
## 11 6.0 0.9000 -1.930 0.978 35.80 18.10 346
## 12 6.5 0.9050 -1.920 0.982 27.10 12.50 294
## 13 7.0 0.9090 -1.900 0.985 20.80 8.71 253
## 14 7.5 0.9080 -1.890 0.984 16.20 6.22 218
## 15 8.0 0.9030 -1.880 0.983 12.80 4.48 190
## 16 8.5 0.9030 -1.860 0.985 10.20 3.25 166
## 17 9.0 0.8990 -1.860 0.984 8.28 2.39 146
## 18 9.5 0.8940 -1.840 0.983 6.75 1.77 129
## 19 10.0 0.8940 -1.820 0.984 5.56 1.32 114
Plot the results.
sizeGrWindow(9, 5)
par(mfrow = c(1,2));
cex1 = 0.9;
# Scale-free topology fit index as a function of the soft-thresholding power
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
xlab="Soft Threshold (power)",
ylab="Scale Free Topology Model Fit,signed R^2",
type="n",
main = paste("Scale independence")
)
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
labels=powers,
cex=cex1,
col="red"
)
# This line corresponds to using an R^2 cut-off
abline(h=0.8,col="red")
# Mean connectivity as a function of the soft-thresholding power
plot(sft$fitIndices[,1], sft$fitIndices[,5],
xlab="Soft Threshold (power)",
ylab="Mean Connectivity",
type="n",
main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5],
labels=powers,
cex=cex1,
col="red"
)
Select best soft power threshold using the above plot. I used a scale-free topology fit index R^2 of 0.831. This lowest recommended R^2 by Langfelder and Horvath is 0.8. I chose 0.83 because we want to use the smallest soft thresholding power that maximizes with model fit. It appears that our soft thresholding power is 3.5 because it is the lowest power before the R^2=0.8 threshold that maximizes with model fit.
Set softPower
softPower <- 3.5
The next few steps may need to be executed on a supercomputer as our dataset is too large for most standard laptops to handle.
# Save data (on local machine)
#save(datExpr, file = "2a-WGCNA.RData")
# Load data (on server)
#adjTOM <- load(file="2a-WGCNA.RData")
May have to execute an the HPC server
Co-expression similarity and adjacency, using the soft thresholding
power saved in softPower and translate the adjacency into
topological overlap matrix to calculate the corresponding dissimilarity.
I will use a signed network. https://peterlangfelder.com/2018/11/25/__trashed/
adjacency <- adjacency(datExpr, power=softPower,type="signed") # Calculate adjacency
TOM <- TOMsimilarity(adjacency, TOMType = "signed") # Translate adjacency into topological overlap matrix
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.
dissTOM <- 1-TOM # Calculate dissimilarity in TOM
Form distance matrix.
geneTree <- flashClust(as.dist(dissTOM), method="average")
We will now plot a dendrogram of genes. Each leaf corresponds to a gene, branches grouping together densely are interconnected, highly co-expressed genes.
plot(geneTree, xlab="", sub="", main="Gene Clustering on TOM-based dissimilarity", labels= FALSE, hang=0.04)
Module identification is essentially cutting the branches off the
tree in the dendrogram above. We like large modules, so we set the
minimum module size (minModuleSize at the
top of script) relatively high. 30 is a good starting point. Module 0 is
reserved for unassigned genes. The are other modules will be listed
largest to smallest.
dynamicMods <- cutreeDynamic(dendro = geneTree,
distM = dissTOM,
deepSplit = 2,
pamRespectsDendro = FALSE,
minClusterSize = minModuleSize
)
## ..cutHeight not given, setting it to 0.891 ===> 99% of the (truncated) height range in dendro.
## ..done.
table(dynamicMods) # List modules and respective sizes
## dynamicMods
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
## 4315 3886 1409 869 673 529 449 365 338 333 311 237 214 203 175 170
## 17 18 19 20 21 22 23 24 25
## 146 142 139 136 134 98 74 66 56
Save results and reload on local machine if needed.
# Save data (on server)
#save(dynamicMods, geneTree, file = "2a-WGCNA.dyMod_geneTree.RData")
# Load data (on local machine)
#dyMod_geneTree <- load(file = "2a-WGCNA.dyMod_geneTree.RData")
Plot the module assignment under the gene dendrogram.
dynamicColors <- labels2colors(dynamicMods) # Convert numeric labels into colors
table(dynamicColors)
## dynamicColors
## black blue brown cyan darkgreen
## 449 3886 1409 203 98
## darkgrey darkred darkturquoise green greenyellow
## 66 134 74 673 311
## grey60 lightcyan lightgreen lightyellow magenta
## 146 170 142 139 338
## midnightblue orange pink purple red
## 175 56 365 333 529
## royalblue salmon tan turquoise yellow
## 136 214 237 4315 869
plotDendroAndColors(geneTree, dynamicColors, "Dynamic Tree Cut", dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05, main = "Gene dendrogram and module colors")
Plot module similarity based on eigengene value.
# Calculate eigengenes
MEList = moduleEigengenes(datExpr, colors = dynamicColors, softPower = 4.5)
MEs = MEList$eigengenes
# Calculate dissimilarity of module eigengenes
MEDiss = 1-cor(MEs)
# Cluster again and plot the results
METree = flashClust(as.dist(MEDiss), method = "average")
plot(METree, main = "Clustering of module eigengenes", xlab = "", sub = "")
ncol(MEs)
## [1] 25
Prepare trait data. Data has to be numeric, so I will substitute time_points and type for numeric values.
treatmentinfo$num <- c("1")
allTraits <- as.data.frame(pivot_wider(treatmentinfo, names_from = !!rlang::sym(group), values_from = num, id_cols = sample_id))
allTraits[is.na(allTraits)] <- c("0")
datTraits <- allTraits %>% column_to_rownames("sample_id")
Define numbers of genes and samples.
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)
Correlations of traits with eigengenes.
moduleTraitCor = cor(MEs, datTraits, use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples);
Colors=sub("ME","",names(MEs))
moduleTraitTree = hclust(dist(t(moduleTraitCor)), method = "average");
plot(moduleTraitTree, main = paste("'",group,"'"," clustering based on module-trait correlation", sep=''), sub="", xlab="", cex.lab = 1.5, cex.axis = 1.5, cex.main = 2)
Correlations of genes with eigengenes.
moduleGeneCor <- cor(MEs,datExpr)
moduleGenePvalue <- corPvalueStudent(moduleGeneCor, nSamples)
Represent module trait correlations as a heatmap.
textMatrix <- paste(signif(moduleTraitCor, 2), "\n(",signif(moduleTraitPvalue, 1), ")", sep = "")
dim(textMatrix) <- dim(moduleTraitCor)
head(textMatrix)
## [,1] [,2] [,3] [,4]
## [1,] "-0.048\n(0.6)" "0.14\n(0.1)" "0.047\n(0.6)" "0.072\n(0.4)"
## [2,] "-0.13\n(0.2)" "0.066\n(0.5)" "0.15\n(0.1)" "0.13\n(0.2)"
## [3,] "-0.0056\n(1)" "0.0022\n(1)" "-0.15\n(0.1)" "-0.25\n(0.006)"
## [4,] "-0.069\n(0.5)" "0.18\n(0.06)" "0.096\n(0.3)" "0.018\n(0.8)"
## [5,] "0.032\n(0.7)" "0.11\n(0.2)" "-0.66\n(4e-16)" "-0.14\n(0.1)"
## [6,] "0.22\n(0.02)" "0.16\n(0.09)" "-0.3\n(9e-04)" "-0.43\n(1e-06)"
## [,5] [,6] [,7] [,8]
## [1,] "-0.21\n(0.02)" "-0.08\n(0.4)" "0.091\n(0.3)" "-0.12\n(0.2)"
## [2,] "-0.041\n(0.7)" "-0.092\n(0.3)" "0.033\n(0.7)" "-0.33\n(3e-04)"
## [3,] "0.13\n(0.2)" "0.55\n(2e-10)" "-0.045\n(0.6)" "-0.32\n(4e-04)"
## [4,] "-0.23\n(0.01)" "-0.039\n(0.7)" "0.086\n(0.4)" "-0.15\n(0.1)"
## [5,] "0.075\n(0.4)" "0.69\n(7e-18)" "-0.039\n(0.7)" "-0.013\n(0.9)"
## [6,] "0.12\n(0.2)" "0.56\n(4e-11)" "-0.065\n(0.5)" "-0.081\n(0.4)"
## [,9]
## [1,] "-0.047\n(0.6)"
## [2,] "0.087\n(0.4)"
## [3,] "0.14\n(0.1)"
## [4,] "-0.081\n(0.4)"
## [5,] "0.081\n(0.4)"
## [6,] "0.043\n(0.6)"
labeledHeatmap(Matrix = moduleTraitCor,
xLabels = names(datTraits),
yLabels = names(MEs),
ySymbols = names(MEs),
cex.lab.y= 0.55,
cex.lab.x= 0.55,
colors = blueWhiteRed(50),
textMatrix = textMatrix,
setStdMargins = TRUE,
cex.text = 0.25,
textAdj = 0,
zlim = c(-1,1),
main = paste("Module-trait relationships")
)
Represent module trait correlations as a complexHeatmap
# Create list of p-values for eigengene correlation with specific life stages
heatmappval <- signif(moduleTraitPvalue, 1)
# Make list of heatmap row colors
htmap.colors <- names(MEs)
htmap.colors <- gsub("ME", "", htmap.colors)
ht=Heatmap(moduleTraitCor,
name = "Eigengene",
column_title = "Module-Trait Eigengene Correlation",
col = blueWhiteRed(50),
row_names_side = "left", row_dend_side = "left",
width = unit(5, "in"), height = unit(8.5, "in"),
column_order = 1:ncol(moduleTraitCor),
column_dend_reorder = TRUE,
cluster_columns = hclust(dist(t(moduleTraitCor)), method = "average"),
column_dend_height = unit(0.5, "in"),
cluster_rows = METree, row_gap = unit(2.5, "mm"), border = TRUE,
cell_fun = function(j, i, x, y, w, h, col) {
if(heatmappval[i, j] <= 0.05) {
grid.text(sprintf("%s", heatmappval[i, j]), x, y, gp = gpar(fontsize = 8, fontface = "bold"))
} else {
grid.text(sprintf("%s", heatmappval[i, j]), x, y, gp = gpar(fontsize = 8, fontface = "plain"))
}},
column_names_gp = gpar(fontsize = 10),
row_names_gp = gpar(fontsize = 10, alpha = 0.75, border = TRUE, fill = htmap.colors)
)
draw(ht)
We quantify associations of individual genes with life stage by
defining Gene Significance GS as the absolute value of the correlation
between the gene and the trait in group. For each module,
we also define a quantitative measure of module membership MM as the
correlation of the module eigengene and the gene expression profile.
Define variable weight containing the weight column of datTrait.
trait <- as.data.frame(as.numeric(as.factor(treatmentinfo[[group]])))
names(trait) = "group"
dim(trait)
## [1] 117 1
Colors of the modules.
modNames <- substring(names(MEs), 3)
geneModuleMembership <- as.data.frame(cor(datExpr, MEs, use = "p"))
MMPvalue <- as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples))
names(geneModuleMembership) <- paste("MM", modNames, sep="")
names(MMPvalue) <- paste("p.MM", modNames, sep="")
geneTraitSignificance <- as.data.frame(cor(datExpr, trait, use = "p"))
GSPvalue <- as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples))
names(geneTraitSignificance) <- paste("GS.", names(trait), sep="")
names(GSPvalue) <- paste("p.GS.", names(trait), sep="")
Create the starting data frame.
geneInfo0 = data.frame(
gene_id = names(datExpr),
moduleColor = dynamicColors,
geneTraitSignificance,
GSPvalue
)
Order modules by their significance for trait.
modOrder = order(-abs(cor(MEs, trait, use = "p")))
Add module membership information in the chosen order.
for (mod in 1:ncol(geneModuleMembership))
{
oldNames = names(geneInfo0)
geneInfo0 = data.frame(geneInfo0,
geneModuleMembership[, modOrder[mod]],
MMPvalue[, modOrder[mod]]
)
names(geneInfo0) = c(oldNames,
paste("MM.", modNames[modOrder[mod]], sep=""),
paste("p.MM.", modNames[modOrder[mod]], sep="")
)
}
Order the genes in the geneInfo variable first by module color, then by geneTraitSignificance.
geneOrder = order(geneInfo0$moduleColor, -abs(geneInfo0[[paste("p.GS.", names(trait), sep="")]]));
geneInfo = geneInfo0[geneOrder, ]
head(geneInfo)
## gene_id
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 Pocillopora_acuta_KBHIv2___RNAseq.g880.t1
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t Pocillopora_acuta_KBHIv2___RNAseq.9444_t
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1
## moduleColor GS.group
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 black -2.068701e-05
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 black -6.102154e-04
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 black 1.806745e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 black 2.002178e-03
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t black -2.866394e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 black 3.701016e-03
## p.GS.group MM.green p.MM.green
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 0.9998234 0.143501397 0.12269986
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 0.9947902 -0.011848636 0.89910582
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 0.9845754 0.155248958 0.09464378
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 0.9829071 -0.080992897 0.38534321
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t 0.9755312 -0.010973204 0.90652525
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 0.9684097 -0.004831206 0.95877036
## MM.yellow p.MM.yellow
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.34892453 0.0001154571
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 -0.30995655 0.0006717571
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 -0.04759341 0.6103575539
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.08257903 0.3760748770
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t -0.03230084 0.7295486759
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.03111254 0.7391335203
## MM.magenta p.MM.magenta
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.06055089 0.516647736
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 -0.13133985 0.158089258
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 0.28681871 0.001717039
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.07031045 0.451274908
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t -0.19967788 0.030893385
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.02746369 0.768810515
## MM.lightcyan p.MM.lightcyan
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.09183281 0.324751101
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 -0.04236013 0.650203857
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 -0.25115627 0.006306037
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.09893561 0.288567834
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t 0.11024440 0.236697001
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 0.08171378 0.381113873
## MM.darkgrey p.MM.darkgrey
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.39969467 8.026113e-06
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 -0.25671100 5.207297e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 -0.01807584 8.466170e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.09963858 2.851386e-01
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t -0.02499447 7.890876e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.09713966 2.974529e-01
## MM.brown p.MM.brown
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.48062039 4.150939e-08
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 -0.50966309 4.389294e-09
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 -0.11151680 2.312966e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.17399406 6.063274e-02
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t -0.02069013 8.247663e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.12194364 1.902791e-01
## MM.darkgreen p.MM.darkgreen
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.20760835 0.02470093
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 -0.09794059 0.29346839
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 0.07561684 0.41776749
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 0.04357657 0.64084788
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t -0.01772474 0.84956079
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.01279809 0.89106891
## MM.royalblue p.MM.royalblue
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.07197915 0.440579294
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 -0.13942809 0.133799310
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 -0.29455654 0.001265272
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.09559146 0.305255539
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t 0.15729986 0.090312512
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 0.05430192 0.560911559
## MM.darkturquoise p.MM.darkturquoise
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.33623611 2.102125e-04
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 -0.37199972 3.630637e-05
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 0.01426237 8.786959e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.30247252 9.177075e-04
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t -0.24991345 6.578395e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.13125934 1.583466e-01
## MM.darkred p.MM.darkred
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.38741657 1.593677e-05
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 -0.30741407 7.475506e-04
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 0.02616069 7.794922e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.06924199 4.581987e-01
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t -0.21558140 1.957728e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.06917848 4.586121e-01
## MM.midnightblue p.MM.midnightblue
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.57279046 1.498926e-11
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 -0.48048974 4.191157e-08
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 -0.04847061 6.037845e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.36559497 5.050166e-05
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t -0.29965632 1.029850e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.23009610 1.256980e-02
## MM.grey60 p.MM.grey60
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.49092055 1.916181e-08
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 -0.57232177 1.570571e-11
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 -0.05451036 5.594064e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.27574594 2.619439e-03
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t -0.12113522 1.932563e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.17041330 6.621404e-02
## MM.lightgreen p.MM.lightgreen
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.30367702 0.0008732472
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 -0.32251816 0.0003905815
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 0.18435207 0.0466183089
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.12317429 0.1858110291
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t -0.26891663 0.0033705434
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.09830956 0.2916447722
## MM.purple p.MM.purple
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 0.6406916 7.281856e-15
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 0.3125822 6.009407e-04
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 0.2499209 6.576734e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 0.2765440 2.542343e-03
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t 0.1864207 4.416974e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 0.3168960 4.993218e-04
## MM.greenyellow p.MM.greenyellow
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.6903599 7.255620e-18
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 -0.6298564 2.791317e-14
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 -0.1537151 9.798913e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.2641491 4.004189e-03
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t -0.2392598 9.375404e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.2517965 6.169689e-03
## MM.lightyellow p.MM.lightyellow
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 0.6421866 6.024165e-15
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 0.3666212 4.792289e-05
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 0.2358117 1.048198e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 0.2110167 2.238515e-02
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t 0.1574486 9.000466e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 0.3705038 3.923999e-05
## MM.pink p.MM.pink MM.tan
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.20804558 0.02439285 -0.5528863
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 -0.18357048 0.04757232 -0.4514506
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 0.20585518 0.02596985 -0.1094448
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 0.05946854 0.52418457 -0.1786946
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t -0.17070584 0.06574293 -0.1504302
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.02806364 0.76390665 -0.1746129
## p.MM.tan MM.red p.MM.red
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 1.022145e-10 0.6079753 3.607703e-13
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 3.244156e-07 0.4803911 4.221761e-08
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 2.401354e-01 0.1739004 6.077368e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 5.389617e-02 0.4355578 9.199729e-07
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t 1.054672e-01 0.1539702 9.742643e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 5.970825e-02 0.2706109 3.168059e-03
## MM.turquoise p.MM.turquoise
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 0.3097140 0.0006786718
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 0.3345102 0.0002276165
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 0.2110974 0.0223326587
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 0.2348850 0.0107982560
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t 0.1835526 0.0475943687
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 0.1292410 0.1649004158
## MM.orange p.MM.orange MM.blue
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.18192643 0.049631983 0.07822197
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 -0.25421504 0.005677917 -0.17397229
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 0.16746756 0.071112581 0.12847846
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 0.18236747 0.049072337 0.06880252
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t -0.02848361 0.760479455 0.07526046
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.14915004 0.108499564 0.09219025
## p.MM.blue MM.black p.MM.black
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 0.40186157 0.7200100 5.803009e-20
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 0.06066546 0.5701121 1.955352e-11
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 0.16742817 0.2922912 1.384786e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 0.46106352 0.5640874 3.524952e-11
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t 0.41997154 0.2995879 1.032722e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 0.32286342 0.3053941 8.132671e-04
## MM.cyan p.MM.cyan MM.salmon
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 0.4415957 6.230095e-07 -0.60650930
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 0.3530400 9.453814e-05 -0.68246587
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 0.2194371 1.744741e-02 -0.08393264
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 0.5526485 1.045064e-10 -0.36854466
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t 0.2462552 7.441642e-03 -0.16963764
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 0.1505055 1.052907e-01 -0.11927007
## p.MM.salmon
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 4.252695e-13
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 2.383775e-17
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 3.682737e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 4.341841e-05
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t 6.747649e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 2.002537e-01
Show and save geneInfo as a TSV.
dim(geneInfo)
## [1] 15467 54
head(geneInfo)
## gene_id
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 Pocillopora_acuta_KBHIv2___RNAseq.g880.t1
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t Pocillopora_acuta_KBHIv2___RNAseq.9444_t
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1
## moduleColor GS.group
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 black -2.068701e-05
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 black -6.102154e-04
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 black 1.806745e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 black 2.002178e-03
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t black -2.866394e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 black 3.701016e-03
## p.GS.group MM.green p.MM.green
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 0.9998234 0.143501397 0.12269986
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 0.9947902 -0.011848636 0.89910582
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 0.9845754 0.155248958 0.09464378
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 0.9829071 -0.080992897 0.38534321
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t 0.9755312 -0.010973204 0.90652525
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 0.9684097 -0.004831206 0.95877036
## MM.yellow p.MM.yellow
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.34892453 0.0001154571
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 -0.30995655 0.0006717571
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 -0.04759341 0.6103575539
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.08257903 0.3760748770
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t -0.03230084 0.7295486759
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.03111254 0.7391335203
## MM.magenta p.MM.magenta
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.06055089 0.516647736
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 -0.13133985 0.158089258
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 0.28681871 0.001717039
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.07031045 0.451274908
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t -0.19967788 0.030893385
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.02746369 0.768810515
## MM.lightcyan p.MM.lightcyan
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.09183281 0.324751101
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 -0.04236013 0.650203857
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 -0.25115627 0.006306037
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.09893561 0.288567834
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t 0.11024440 0.236697001
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 0.08171378 0.381113873
## MM.darkgrey p.MM.darkgrey
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.39969467 8.026113e-06
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 -0.25671100 5.207297e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 -0.01807584 8.466170e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.09963858 2.851386e-01
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t -0.02499447 7.890876e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.09713966 2.974529e-01
## MM.brown p.MM.brown
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.48062039 4.150939e-08
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 -0.50966309 4.389294e-09
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 -0.11151680 2.312966e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.17399406 6.063274e-02
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t -0.02069013 8.247663e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.12194364 1.902791e-01
## MM.darkgreen p.MM.darkgreen
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.20760835 0.02470093
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 -0.09794059 0.29346839
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 0.07561684 0.41776749
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 0.04357657 0.64084788
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t -0.01772474 0.84956079
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.01279809 0.89106891
## MM.royalblue p.MM.royalblue
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.07197915 0.440579294
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 -0.13942809 0.133799310
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 -0.29455654 0.001265272
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.09559146 0.305255539
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t 0.15729986 0.090312512
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 0.05430192 0.560911559
## MM.darkturquoise p.MM.darkturquoise
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.33623611 2.102125e-04
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 -0.37199972 3.630637e-05
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 0.01426237 8.786959e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.30247252 9.177075e-04
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t -0.24991345 6.578395e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.13125934 1.583466e-01
## MM.darkred p.MM.darkred
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.38741657 1.593677e-05
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 -0.30741407 7.475506e-04
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 0.02616069 7.794922e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.06924199 4.581987e-01
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t -0.21558140 1.957728e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.06917848 4.586121e-01
## MM.midnightblue p.MM.midnightblue
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.57279046 1.498926e-11
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 -0.48048974 4.191157e-08
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 -0.04847061 6.037845e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.36559497 5.050166e-05
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t -0.29965632 1.029850e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.23009610 1.256980e-02
## MM.grey60 p.MM.grey60
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.49092055 1.916181e-08
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 -0.57232177 1.570571e-11
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 -0.05451036 5.594064e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.27574594 2.619439e-03
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t -0.12113522 1.932563e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.17041330 6.621404e-02
## MM.lightgreen p.MM.lightgreen
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.30367702 0.0008732472
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 -0.32251816 0.0003905815
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 0.18435207 0.0466183089
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.12317429 0.1858110291
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t -0.26891663 0.0033705434
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.09830956 0.2916447722
## MM.purple p.MM.purple
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 0.6406916 7.281856e-15
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 0.3125822 6.009407e-04
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 0.2499209 6.576734e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 0.2765440 2.542343e-03
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t 0.1864207 4.416974e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 0.3168960 4.993218e-04
## MM.greenyellow p.MM.greenyellow
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.6903599 7.255620e-18
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 -0.6298564 2.791317e-14
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 -0.1537151 9.798913e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.2641491 4.004189e-03
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t -0.2392598 9.375404e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.2517965 6.169689e-03
## MM.lightyellow p.MM.lightyellow
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 0.6421866 6.024165e-15
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 0.3666212 4.792289e-05
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 0.2358117 1.048198e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 0.2110167 2.238515e-02
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t 0.1574486 9.000466e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 0.3705038 3.923999e-05
## MM.pink p.MM.pink MM.tan
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.20804558 0.02439285 -0.5528863
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 -0.18357048 0.04757232 -0.4514506
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 0.20585518 0.02596985 -0.1094448
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 0.05946854 0.52418457 -0.1786946
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t -0.17070584 0.06574293 -0.1504302
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.02806364 0.76390665 -0.1746129
## p.MM.tan MM.red p.MM.red
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 1.022145e-10 0.6079753 3.607703e-13
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 3.244156e-07 0.4803911 4.221761e-08
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 2.401354e-01 0.1739004 6.077368e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 5.389617e-02 0.4355578 9.199729e-07
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t 1.054672e-01 0.1539702 9.742643e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 5.970825e-02 0.2706109 3.168059e-03
## MM.turquoise p.MM.turquoise
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 0.3097140 0.0006786718
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 0.3345102 0.0002276165
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 0.2110974 0.0223326587
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 0.2348850 0.0107982560
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t 0.1835526 0.0475943687
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 0.1292410 0.1649004158
## MM.orange p.MM.orange MM.blue
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.18192643 0.049631983 0.07822197
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 -0.25421504 0.005677917 -0.17397229
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 0.16746756 0.071112581 0.12847846
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 0.18236747 0.049072337 0.06880252
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t -0.02848361 0.760479455 0.07526046
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.14915004 0.108499564 0.09219025
## p.MM.blue MM.black p.MM.black
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 0.40186157 0.7200100 5.803009e-20
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 0.06066546 0.5701121 1.955352e-11
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 0.16742817 0.2922912 1.384786e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 0.46106352 0.5640874 3.524952e-11
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t 0.41997154 0.2995879 1.032722e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 0.32286342 0.3053941 8.132671e-04
## MM.cyan p.MM.cyan MM.salmon
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 0.4415957 6.230095e-07 -0.60650930
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 0.3530400 9.453814e-05 -0.68246587
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 0.2194371 1.744741e-02 -0.08393264
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 0.5526485 1.045064e-10 -0.36854466
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t 0.2462552 7.441642e-03 -0.16963764
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 0.1505055 1.052907e-01 -0.11927007
## p.MM.salmon
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 4.252695e-13
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1 2.383775e-17
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1 3.682737e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 4.341841e-05
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t 6.747649e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 2.002537e-01
write.table(geneInfo, file=WGCNA_results.file, sep='\t', quote=F, row.names=F)
write.table(treatmentinfo, file=WGCNA_samples.file, sep='\t', quote=F, row.names=F)